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Preface 


This  thesis  began  as  an  attempt  to  duplicate  the  suspect  results  of 
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implicit  method  of  Peaceman  and  Rachford.  A  literature  review  has  shown 
that  this  subject  has  not  been  addressed  previously  in  the  literature. 
Chiver's  thesis  covers  much  of  the  textbook  material  on  finite 
differences  and  may  be  of  interest  to  some  readers  for  a  review  of  some 
basics. 
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for  sponsoring  this  study  and  to  Sheryl  Michel  for  her  excellent  typing. 
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Notation  and  Abbreviations 


Ax  3  Ay  »  h  3  mesh  spacing,  the  increment  of  space  between  nodal  points, 
t  =  time  (seconds) 

At  =  time  increment 

a  =  thermal  diffusivity ,  here  always  =  1.0 
r  =  aAt/AX*  =  aAt/Ay*  =  aAt/h2 

T(i,j,n)  =  Temperature  (degrees)  at  the  ith,  jth  nodal  point  at  the  nth 
time  step,  3  T( iAx ,jAy ,nAt)  =  T  (x,y,t) 

ADI  *  Peaceman-Rachford  alternating  direction  implicit  method 

CN  =  Crank-Nicolson  implicit  method 

IM  =  Fully  implicit  method 

EX  =  Fully  explicit  method 

_A  3  matrix  notation  for  the  matrix  A 

b  =  vector  notation  for  the  vector  b 


Abstract 
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^yThe  two  dimensional  heat  diffusion  equation  with  Dirichlet  boundary 
conditions  was  solved  using  the  fully  explicit,  fully  implicit, 
Crank-Nicolson  implicit,  and  Peaceman-Rachford  alternating  direction 

r 

implicit  (ADI)  methods.  Comparisons  of  accuracy  and  time  requirements 
were  made. 

The  possibility  that  the  ADI  method  has  stable  oscillatory 
solutions  with  large  time  steps  was  investigated. 

Results  of  computations  revealed  that  the  ADI  has  stable 
oscillations  for  large  time  steps,  in  some  cases  producing  large  enough 
errors  to  render  the  solution  unusable.  Time  steps  greater  than  twice 
the  square  of  the  mesh  spacing  divided  by  the  thermal  diffusivity  must 
be  used  with  care. 

For  small  time  steps,  the  Crank-Nicolson  and  ADI  methods  were  the 
most  accurate,  and  the  ADI  was  the  fastest  method.  The  fully  implicit 
method  was  the  most  accurate  at  large  time  steps,  but  the  ADI,  with  a 
smaller  time  step  to  reduce  the  oscillatory  error,  was  still  the  fastest 


method  to  reach  a  solution  with  the  desired  degree  of  accuracy. 
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I.  Introduction 


Background 

The  heat  diffusion  equation  is  encountered  in  many  areas  of 
engineering,  including  nuclear  engineering.  Since  only  the  simplest 
problems  have  analytical  solutions,  numerical  techniques  utilizing 
computers  are  commonly  used  to  solve  the  practical  problems  encountered 
in  engineering.  Finite  difference  methods  are  the  most  basic  and 
widespread  techniques  for  obtaining  solutions  to  the  heat  diffusion 
equation. 

Although  numerous  schemes  appear  in  the  literature,  the  standard 
finite  difference  methods  used  in  solving  the  heat  diffusion  equation 
remain  the  fully  explicit  (EX),  fully  implicit  (IM),  Crank-Nicol son 
implicit  (CN),  and  the  Peaceman-Rachford  alternating  direction  implicit 
(ADI)  methods.  The  employment  of  these  methods  requires  an 
understanding  of  their  behavior  under  different  conditions,  since  large 
errors  may  sometimes  occur. 


The  Two  Dimensional  Heat  Diffusion  Equation 

The  heat  diffusion  equation  in  two  dimensions  is  a  parabolic 
partial  differential  equation  of  the  form 


3T(x,y,t) 

at 


«[■ 


n(x,y,t) 

3x2 


3  T(x,y,t) 

1/ 


(1) 


where  T  is  the  temperature  and  o  is  the  thermal  diffusivity.  This 
equation  describes  the  diffusion  of  heat  in  a  solid  in  which  the  heat 
diffusion  is  independent  of  the  z  direction  (4:12). 


To  solve  a  problem,  the  initial  and  boundary  conditions  must  be 
specified.  The  three  standard  types  of  boundary  conditions  are  the 
Dirichlet,  which  specifies  the  value  of  T  at  the  boundaries,  the 
Neumann,  which  specifies  the  value  of  the  derivative  of  T  at  the 
boundaries,  and  the  mixed,  which  combines  the  Dirichlet  and  Neumann 
(4:21). 

Problem  Statement 

The  problem  investigated  in  this  thesis  is  to  compare  the  accuracy 
and  time  required  for  the  fully  explicit,  fully  implicit,  Crank- 
Nicolson,  and  Peaceman-Rachford  alternating  direction  implicit  methods 
of  solution  of  the  two  dimensional  heat  diffusion  equation.  Further, 
the  possibility  is  investigated  that  the  ADI  method  has  stable 
oscillations  under  certain  conditions,  similar  to  the  stable 
oscillations  known  to  exist  in  the  Crank-Nicolson  method  (7:284;8:57) . 
Finally,  it  should  be  noted  that  this  thesis  is  a  follow-on  to  a 
previous  thesis  by  Chivers.  His  results,  that  the  ADI  and  CN  methods 
had  unstable  oscillations  (3:35),  do  not  agree  with  theory  and  his 
error,  if  any,  needed  to  be  found. 

Scope 

This  study  was  restricted  to  the  two  dimensional  heat  diffusion 
equation  over  a  unit  square.  The  mesh  was  set  up  sucn  that  Ax  =  Ay  =  h, 
only  Dirichlet  boundary  conditions  were  used,  and  the  thermal 
diffusivity  was  1.0. 


General  Approach 

Initially,  a  previous  thesis  by  Chivers  (3)  was  reviewed  to  try  to 
find  any  errors  in  his  approach,  and  a  literature  review  was  conducted. 
After  these,  the  coefficient  method  was  used  to  determine  the  stability 
of  the  ADI  method  in  two  dimensions.  Computer  programs  for  the  four 
finite  difference  methods  were  then  written  to  solve  the  two  dimensional 
heat  diffusion  equation  with  Dirichlet  boundary  conditions. 

Two  problems  with  known  analytical  solutions,  allowing  comparison, 
were  studied.  The  first  problem  solved  was  the  same  one  Chivers  used, 
since  it  was  desired  to  either  duplicate  his  results  or  determine  his 
error.  A  second  problem  was  then  selected  to  give  more  insight  into  the 
behavior  of  the  finite  difference  methods  with  a  practical  engineering 
problem.  This  problem  was  solved,  after  which  the  oscillations  of  the 
ADI  method  were  investigated  in  more  detail. 

Sequence  of  Presentation 

Chapter  2  presents  general  finite  difference  equations  along  with 
the  specific  formulations  of  the  four  methods.  Systems  of  equations, 
resulting  the  replacement  of  the  heat  equation  with  finite 
differences,  a?\.  o  discussed. 

An  analysis  of  stability  of  the  four  methods,  using  the  coefficient 
method,  is  presented  in  Chapter  3.  Also  discussed  in  this  chapter  are 
the  expected  characteristics  of  the  finite  difference  methods,  such  as 
accuracy  and  time  required. 

Chapter  4  presents  results  of  the  computer  solutions  to  the  heat 
diffusion  equation  using  the  four  different  methods.  Conclusions  and 
recommendations  are  given  in  Chapter  5. 
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is  often  the  slowest  method  (9:27).  The  ADI  may  be  fastest  as  a  result 
of  the  tridiagonal  form  of  its  matrix  (5:111),  and  tne  more  grid  points 
involved,  the  more  advantage  it  has  over  the  CN  and  IM  methods,  as 
discussed  in  Chapter  2.  Whether  the  IM  method  has  an  advantage  over  the 
ADI  for  larger  At  may  depend  on  th  size  of  the  error  created  by  any 
oscillations  of  the  ADI  method.  If  forced  to  use  a  smaller  At  due  to 
oscillations,  the  method  could  conceivably  take  longer  than  the  IM 
method.  This  is  investigated  in  the  next  chapter. 
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where 

I  3  the  number  of  interior  nodal  points  in  x 
J  *  the  number  of  interior  nodal  points  in  y 

Ta(i,j,n)  3  the  analytical  solution  at  the  ith,  j:h  grid  point  at  the 
nth  time  step 

T(i.j.n)  3  the  computed  solution  at  the  ith,  jth  grid  point  at  the 
nth  time  step 

This  uses  the  standard  definition  of  root  mean  squared  average  (1:506) 
and  varies  slightly  from  the  one  used  by  Towler  and  Yang  (  14:1023)  in 
that  IJ  is  under  the  square  root.  Note  that  since  Oirichlet  boundary 
conditions  are  used,  the  boundary  nodal  points  are  not  included  in  the 
average  (14:1023). 

Besides  using  an  average  error,  the  maximum  error  at  any  nodal 
point  is  a  useful  tool,  since  a  small  average  does  not  insure  accuracy 
at  all  the  points.  A  determination  of  the  maximum  error  shows  whether 
the  solution  has  accuracy  over  the  entire  grid,  which  is  important  for 
most  engineering  applications. 

Computer  T ime 

For  a  set  number  of  time  steps,  the  EX  method  should  be  fastest, 
and  the  CN  and  IM  slowest.  This  is  because  the  EX  solves  for  the 
temperature  directly  at  each  nodal  point.  The  ADI,  as  previously  noted, 
solves  a  system  of  equations  with  a  tridiagonal  matrix,  so  that  it  is 
faster  per  time  step  than  the  IM  and  CN  methods. 

For  actually  solving  a  problem,  the  stability  restriction  on  the  EX 
method  forces  many  times  steps  with  a  small  At  so  that  the  EX  method 
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Accuracy  of  Methods 

The  accuracy  of  these  methods  due  to  truncation  error  of  the  Taylor 
series  approximations  are  of  the  order  At*  +  h*  for  the  AO  I  and  CN 
methods,  and  of  the  order  At  +  h*  for  the  EX  and  IM  methods 
(11:212; 13:42 ,90) .  Generally,  the  numerical  solutions  are  more  accurate 
than  these  estimates.  Since  not  all  the  errors  have  the  same  sign,  they 
tend  to  partially  cancel  (13:57). 

Other  errors  are  introduced  besides  truncation  error.  Oscillatory 
solutions  of  the  CN  method  -  and  as  proposed  here  of  the  ADI  -  for  some 
values  of  r  may  add  significant  error.  Oscillatory  behavior  of  the  EX 
method  for  r  >  1/4  results  in  gross  error.  Round  off  error  by  the 
computer  may  increase  when  a  large  number  of  arithmetic  operations  are 
required,  as  when  solving  a  system  of  simultaneous  equations. 

Finally,  it  should  be  noted  it  is  widely  known  that  the  CN  method 
is  more  accurate  than  the  I M  method  at  small  time  steps,  while  the  IM 
method  is  more  accurate  than  the  CN  method  for  large  time  steps 
(9:27-28). 

Measuring  Accuracy 

The  accuracy  of  the  solutions  was  found  by  computing  the  root  mean 
square  (RMS)  average  error  of  all  the  grid  points  and  by  finding  the 
maximum  error  at  any  of  the  grid  points.  These  errors  were  computed 
using  the  analytical  and  computed  solutions  at  each  nodal  point.  This 
follows  the  convention  suggested  by  Towler  and  Yang  (14:1023-1024).  The 
RMS  error  is  computed  as 

J  I  1/2 
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Figure  1.  Stability  Curves  For  One  and  Two  Dimensions 
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This  will  always  be  between  the  values  of  *1  and  -1.  Therefore,  the  CN 
formulation  is  always  stable,  but  the  solution  may  oscillate  when  the 
temperature  ratio  is  less  than  zero,  which  occurs  when 


r  >  1/2 


Again,  this  is  more  restrictive  than  the  one  dimensional  cases  where 
stable  oscillations  may  occur  if  r  >  1  (9:30).  Note  that  any 
oscillations  will  eventually  die  out  since  the  error  is  bounded. 

The  ADI  formula,  Eq  (5),  gives,  for  a  one  nodal  point  grid 


TCi»4»n^lJ 

TU.j.n) 


_  1  -  2r 
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This  is  the  same  stability  criteria  as  for  the  CN  method,  Eq  (14).  Thus 
the  ADI  method  is  always  stable  but  may  have  oscillatory  solutions. 

The  stability  curves  of  all  four  methods  for  the  one  and  two 
dimensional  cases  are  plotted  in  Figure  1.  Although  here,  Ax  »  Ay,  the 
above  analysis  can  easily  be  extended  to  a  grid  system  where  Ax  A  Ay. 
For  example,  for  the  ADI  and  CN  method,  Eqs  (14)  and  (16)  become 
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The  possible  ranges  of  the  ratfo  are 


when 


and 


when 


and 


when 


0  <  T( i , j ,n+l)/T( i ,  j  ,n)  <  1 

r  <  1/4 

-1  <  T(i,j,n+1)/T(i,j,n)  <  0 

1/4  <  r  <  1/2 

T(i,j,n+1)/T(i,j,n)  <  -1 

r  >  1/2 


(7) 

(8) 

(9) 

(10) 

(11) 

(12) 


Note  that  this  is  more  restrictive  than  the  one  dimensional  case,  where 
the  stability  condition  for  T( i , j ,n+l)/T( i , j  ,n)  >  0  is  that  r  <  1/2 
(5:108). 

The  fully  implicit  formula,  Eq  (3),  yields 


T(i,j,n) 

which  is  always  positive 
The  CN  formulation. 


1 

1  +  4r 

and,  therefore,  always 

Eq  (4),  yields 

1  -  2r 

i"  “  77 


(13) 

stable. 


(14) 
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1.  0  <  T(i,j,n+1)/T(i,j,n)  <  1.  Stable  solution,  the  error 
decreases  with  time,  no  oscillations. 

2.  -1  <  T(i,j,n+1)/T(i,j,n)  <  0.  Stable  oscillations,  error 
and  oscillations  decrease  with  time. 

3.  T(i,j,n+1)/T(i,j,n)  <  -1.  Unstable  oscillations,  error  is 
unbounded,  growing  with  time. 

Note  that  the  ratio  is  never  greater  than  1. 

This  behavior  is  easily  understood.  If  the  temperature  ratio  is 
positive,  each  time  step  results  in  a  temperature  of  the  same  sign  but 
with  smaller  magnitude.  Thus  the  temperature  steadily  decreases  with 
time,  and  is  stable. 

If  the  temperature  ratio,  is  negative  but  greater  than  -1,  each  new 
temperature  will  have  the  opposite  sign  of  the  previous  temperature  and 
a  smaller  magnitude.  Thus  the  temperature  will  oscillate  between 
positive  and  negative  values,  but  will  decrease  in  magnitude  as  time 
steps  are  taken,  and  so  it  is  stable. 

If  the  temperature  ratio  is  less  than  -I,  each  new  temperature  will 
have  the  opposite  sign  of  the  preceding  temperature  and  will  be  larger 
in  magnitude.  Thus  the  temperature  will  oscillate  while  its  magnitude 
will  grow  without  bound,  causing  instability. 


Stability  for  Difference  Methods  jn  Two  Dimensions 

In  two  dimensions,  the  explicit  formula,  Eq  (2),  evaluated  for  a 
one  nodal  point  grid,  yields 


T( i , j ,n+l) 
T(i,j,n) 


1  -  4r 


(6) 


III.  Stab i 1 ity  and  Properties 


Stability 

The  concept  of  stability  is  sometimes  misunderstood.  Stability  is 
simply  the  guarantee  that  the  error  does  not  grow  in  time,  that  is,  the 
error  is  bounded  (5:98).  Stability  of  a  finite  difference  formulation 
does  not  guarantee  accuracy  or  non-oscillatory  behavior  of  the  solution, 
although  any  oscillations  will  die  out  eventually  (9:27).  On  the  other 
hand,  a  non-stable  formulation  will  blow  up  as  the  number  of  time 
intervals  increase,  that  is,  the  error  will  grow  with  time  in  an 
unbounded  fashion. 

A  common  and  easily  understood  method  of  considering  stability  is 
the  coefficient  method.  Myers  and  Patankar  both  present  this  method, 
but  only  for  the  one  dimensional  heat  diffusion  equation  (7:281-284; 
9:28-30).  The  method  is  developed  here  for  two  dimensions. 

The  idea  of  the  coefficient  method  is  to  consider  a  system 
consisting  of  only  one  nodal  point,  so  that  the  only  non-zero 
temperatures  are  T(i,j,n)  and  T(i,j,n+1).  Then  the  ratio 


Ciiij.iLL, 

t(i»J,n+l) 


where  C(i,j,n)  is  the  coefficient  of  T(i,j,n),  is  easily  found.  The 
three  types  of  behavior  of  the  solution  for  different  values  of  the 
ratio  T(i,j,n+1)/T(i,j,n)  are  (7:282): 


The  IM  and  CN  methods,  for  the  problems  used  in  this  study,  result 
in  being  a  real,  symmetric,  positive  definite  matrix.  The  most 
efficient  method  for  solving  the  equation  Ab=c  for  this  case  is  the 
direct  Cholesky  square  root  method  (15:105).  This  method  requires  on 
the  order  of  N^/3  operations,  although  it  is  faster  than  this  sounds 
since  the  matrix  -A  has  only  5  non-zero  diagonals  (15:99-100).  Should  a 
problem  result  in  a  non-symmetric  matrix,  a  less  efficient  method,  such 
as  the  Gauss-Seidel ,  would  be  needed.  For  this  reason,  a  problem  with 
irregular  geometry  which  does  not  result  in  a  symmetric,  positive 
definite  matrix  will  take  more  work  to  solve  than  the  problems  used  in 
this  study. 

The  ADI  method  results  in  a  tridiagonal  matrix  A..  This  is  a  great 
advantage,  as  the  equation  Ab=c  can  be  solved  rapidly  using  the  Thomas 
method.  Since  this  method  only  uses  on  thejjrder  of  8N  operations,  it 
is  very  fast  (5:111).  Note  that  the  larger  the  number  of  mesh  points  N, 
the  faster  the  ADI  will  be  relative  to  the  IM  and  CN  due  to  the  order  of 
operations  being  8N  vs  N^/3.  Additionally,  a  tridiagonal  matrix  in  band 
storage  mode  requires  much  less  storage,  which  is  a  significant  factor 
for  a  system  with  a  large  number  of  mesh  points. 


(l+2r)  T( i ,j ,n+l)  -  r[T(i+l,j,n+l)  ♦  T(i-l,j,n+l)] 

»  (l-2r)  T(i,j,n)  ♦  r[T( i ,j+l,n)  ♦  T(i,j-l,n)] 


to  advance  from  time  step  n  to  n+1,  followed  by 


d+2r)  T(i,j,n+2)  -  r[T(i,j+l,n+2)  ♦  T(1,j-l,n+2)] 

»  (l-2r)  T(i ,j ,n+l)  +  r [T( 1 ♦! , j . n+1)  ♦  T( 1 -1 , j . n+1) ] 


to  advance  from  time  step  n+1  to  n+2.  Note  that  At  must  remain 
constant,  and  that  the  solution  is  valid  only  after  both  steps  are  taken 
(6:53; 13:42). 


Systems  of  Equations 

The  EX  method  results  In  one  equation  per  mesh  point  to  solve  for 
the  temperature  at  t  ♦  At  directly.  The  other  methods,  however,  result 
in  systems  of  equations  which  must  be  solved  simultaneously.  For  N 
interior  mesh  points,  this  is. represented  in  matrix  form  as 


where  i  is  a  NxN  matrix,  b  is  a  vector  composed  of  the  new  temperature 
at  t  +  At,  and  c  is  a  vector  composed  from  the  old  temperatures  and 
boundary  conditions.  This  equation  is  then  solved  for  b. 


Fully  Imp! icit  (IM)  Formulation 

The  fully  implicit  method  computes  the  temperature  at  the  new  time 
step  based  on  temperatures  at  neighboring  nodal  points  which  are  also  at 
the  new  time.  A  simple  extension  of  the  one  dimensional  case  (8:58) 
results  in  a  series  of  equations  of  the  form 


(l+4r)  T(i,j,n+1)  -  r[T(i+l,j ,n+l) 
♦  T(i,j-l,n+l)]  *  T(i,j,n) 


♦  T(i-l,j,n+l)  +  T(i ,j+l,n+l) 

(3) 


which  must  be  solved  simultaneously. 

Crank-Nicolson  (CN)  Imp! icit  Formulation 

The  Crank-Nicolson  method  weighs  the  fully  explicit  and  fully 
implicit  equally  and  combines  the  two  methods.  The  result  is  a  series 
of  equations  which  again  must  be  solved  simultaneously.  The  equations 
are  of  the  form  (13:42) 


i 


(2+4r)  T( i , j ,n+l)  -  r[T( i+1 , j ,n+l)  + 

♦  T(i,j-l,n+l)] 

*  (2-4r )  T(i,j,n)  +  r[T(i+l,j,n) 
'♦  T(i ,j-l,n)] 


T( i-l»j »n+l)  ♦  T(i,j+l,n+l) 

+  T(i-l,j,n)  +  T(i ,j+l,n) 


(4) 


Peaceman-Rachford  (ADI)  Formulation 

The  ADI  method  alternates  solving  one  space  direction  explicitly 


and  the  other  implicitly.  The  method  results  in  a  series  of 
simultaneous  equations  to  solve  of  the  form  (5:110;1Q:33) 


Standard  textbooks  in  this  field  cover  all  the  variations  from  which 


finite  differences  may  be  formulated  ( 5 :59 ;7  : 234 ;8  : 25- 28 ) .  The  result  ^ , 

,*.v 

of  using  finite  differences  is  that  a  partial  differential  equation  may  ,';X 

be  replaced  by  a  series  of  finite  difference  equations.  The  solution  of  X.* 

the  finite  difference  equations  on  a  computer  then  gives  the  answer  at  tv 

r.V 

each  nodal  point.  X> 


Heat  Diffusion  Equation  l£  Finite  Difference  Form 

The  four  primary  methods  of  solving  the  heat  diffusion  equation 
with  finite  differences  are  the  fully  explicit,  fully  implicit, 
Crank-Nicolson  implicit,  and  Peaceman-Rachford  alternating  direction 
implicit  methods.  They  are  found  in  numerous  texts  and  journals  in  the 
field  of  heat  conduction.  As  a  reminder,  for  this  study,  h*Ax«Ay,o*l, 
and  r«aAt/hl*At/h* . 

Fully  Exp! icit  (EX)  Formulation 

The  explicit  name  implies  the  form  of  this  solution.  It  computes 
temperature  at  each  nodal  point  in  terms  of  the  old  temperatures  of  the 
previous  time  step.  The  formula  is  (13:41), 


T(i,j,n+1)  *  r[T( i+1 ,j ,n)  +  T(i-l,j,n)  ♦  T(i,j+l,n) 
+  T(i,j-l,n)]  ♦  (l-4r)T(i,j,n) 


(2) 


This  method  is  the  most  straight  forward  of  the  four  methods  as  it  only 
requires  the  solution  of  an  independent  algebraic  equation  for  each 
point  for  every  new  time  step. 


f ' 


r-r 
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II.  Finite  Difference  Methods 


Finite  Difference  Equations 

To  formulate  two  dimensional  finite  difference  equations,  time  and 
space  are  divided  into  discrete  intervals  At,  Ax,  and  Ay,  forming  a 
mesh  or  grid  of  nodal  points.  Instead  of  continuous  values,  there  are 
only  discrete  values  of  variables  and  functions,  so  that  x  *  iAx,  y*jAy, 
and  t*nAt  where  i,j,  and  n  are  integers.  Then  a  function  such  as 
temperature,  T(x,y,t),  is  replaced  with  a  function  with  values  only  at 
mesh  points,  T(x,y,t)  *  T( i Ax , j Ay , nA t ) ,  which  notationally  is 
represented  as  T(i,j,n). 

Once  the  grid  has  been  established,  the  partial  differentials  in  an 
equation  are  replaced  with  finite  differences.  These  are  built  from 
Taylor  series  expansions  of  the  function  evaluated  at  a  point,  for 
example,  adding  the  Taylor  series  expansions 

i  2 

T(x  +  Ax)  *  T(x)  +  ax  T  (x)  +  T"(x)  +  0(ax3)  +  ... 

and 

2 

T(x  -  ax)  *  T(x)  -  ax  T'(x)  +  T"(x)  -  0(ax3)  +  ... 

yields  the  central  difference  to  a  second  derivative, 

T’(x>  *  .  T(x  *  ax)  -  2T(x),  ♦  A»)  +  0(4x2) 

dx^  Ax^ 


Computer 

The  computer  used  was  a  VAX  11/780,  which  is  a  32  bit  machine. 
Double  precision  was  used  to  minimize  round  off  error.  Graphs  were 
drawn  with  a  CDC  CYBER  computer  using  DISSPLA.  All  programs  were 
written  in  FORTRAN  77. 


'  /■  ,*•  / 
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IV.  Results 


This  chapter  presents  the  results  of  numerous  computer  runs  to 
compare  the  behavior  and  characteristics  of  the  different  finite 
difference  methods.  First  the  two  problems  and  their  analytical 
solutions  are  presented.  Then  the  results  for  accuracy,  oscillations, 
and  time  are  presented.  Finally,  a  brief  summation  of  errors  in 
Chiver's  thesis  is  given. 

Problem  Description 

To  compare  the  four  finite  difference  methods,  two  problem’s  with 
known  analytical  solutions  were  solved.  After  getting  results  for  the 
first  problem,  it  became  apparent  that  the  problem  was  inadequate  for 
evaluating  the  usefulness  of  the  methods  for  a  practical  engineering 
problem.  This  was  because  the  solution  rapidly  tends  toward  zero,  and 
the  errors  and  oscillations  were  so  small  that  their  interpretation  was 
difficult.  A  second  problem  which  more  closely  resembles  a  practical 
engineering  problem  was  then  selected  and  solved. 

Problem  #1 

The  first  problem  is 


in  the  region  0  <  x  <  1,  0<y<l,  and  0  <  t,  with  the  initial 
condition 


T(x,y,0)  *  sin(irx)  sin(iry) 

and  boundary  conditions 

T(0,y,t)  *  T(l,y,t)  =  T(x,0,t)  *  T(x,l,t)  *  0 

The  analytical  solution  is 

T(x,y,t)  *  exp(-2ir*t)  sin(*x)  sin(Try) 

which  is  easily  found  using  separation  of  variables  and  verified  by 
substitution.  The  solution  found  by  Chivers  (3:5,55)  is  erroneous  in 
that  the  initial  condition  was  not  properly  applied.  The  problem  is 
similar  to  one  in  Mitchell  (6:66)  except  for  the  boundary  conditions, 
and  the  solution  is  the  same.  For  purposes  of  comparison,  the  solution 
was  computed  for  a  time  of  1.0  seconds. 

Problem  12 

The  second  problem  is 

3T(x,y,t)  B  a2  T(x,y,t)  +  a2  T(x,y,t) 

at  T2  ri 

*  y 
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in  the  region  0  <  x  <  1,  0  <  y  <  1,  and  0  <  t,  with  the  initial 
condition 

T(x,y,0)  -  0 

and  boundary  conditions 

T(O.y.t)  *  T{l,y,t)  -  T(x ,0,t)  -  T(x,l,t)  *  400.0 

This  problem  and  its  analytical  solution  is  virtually  identical  to  the 
one  presented  by  Shih  and  Skladany  (12:410)  except  that  here  it  is  not 
normalized.  The  analytical  solution  is 

T(x ,y,t)  -  400(l-f (x,y,t)) 

f(x,y,t)  =  g(x,t)h(y,t) 

®  4 

g(x,t)  =  z  - —  exp(-n*ir*t)  sin(mrx) 

n=l  nir 

\ 

4 

h(y,t)  *  Z  -  exp(-nlirzt)  sin(niry) 

n-1  nir 

n  — 

The  solution  is  found  by  first  normalizing,  then  using  separation  of 
variables  (2:173;12:41Q) ,  and  can  be  verified  by  substitution.  For 
accuracy  comparisons,  the  solution  was  computed  at  times  of  0.1  and  0.4 
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where 


with 


and 


where 


seconds. 


Figures  2  and  3  present  the  accuracy  results  for  problem  #1,  and 
Figures  4-7  give  the  results  for  problem  #2.  Each  figure  has  two 
graphs,  one  plotting  the  RMS  error  and  one  plotting  the  maximum  error 
for  different  at.  In  Figures  2,  4,  and  6,  the  mesh  spacing  is  h  *  0.1, 
and  in  Figures  3,  5,  and  7,  the  mesh  spacing  is  h  *  0.05,  so  that  the 
effect  of  halving  the  mesh  size  can  be  seen. 

Results 

The  results  of  all  the  cases  in  Figures  2-7  show  that  for  small 
time  steps  the  AOI  and  CN  methods  are  more  accurate  than  the  IM  method. 
This  is  expected  since  they  have  less  truncation  error  (see  Chapter  2). 
At  large  time  steps  the  IM  becomes  more  accurate,  due  to  oscillatory 
solutions  for  the  CN  and  AOI  methods  at  large  time  steps.  The 
oscillations  are  discussed  in  more  detail  in  the  next  section.  Notice 
that  the  CN  and  ADI  solutions  are  very  close  in  accuracy. 

Based  on  truncation  error,  the  EX  method  is  expected  to  be  about  as 
accurate  as  the  IM  method.  As  Figures  2,  6,  and  7  show,  however,  there 
are  regions  where  it  is  the  most  accurate.  This  occurs  with  small  At  , 
which  require  many  times  steps  to  reach  a  solution  at  the  desired  time. 
The  CN,  ADI,  and  IM  methods,  requiring  solutions  of  simultaneous 
equations  for  each  iteration,  accumulate  more  round  off  error,  giving 
the  EX  method  an  overall  accuracy  advantage  in  these  cases.  Note  that 
the  EX  method  gives  unbounded  error  when  r  >  1/4. 
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Figure  3.  Error  vs  Time  Step  Problem  II 
Time«1.0  Ax=Ay*h*0.05 
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TIME  STEP  (seconds) 


Figure  4.  Error  vs  Time  Step  Problem  12 


Time»0.1  Ax»Ayah»0.1 


24 


Figure  7.  Error  Vs  Time  Step  Problem  #2  Time=0.4  Ax=Aysh=0.05 


TIME  STEP  (seconds)  TIME  STEP  (seconds 


For  small  time  steps  (10-^  to  10"^},  the  accuracy  of  the  ADI  and  CN 


methods  remain  about  the  same,  while  the  IM  and  EX  methods  have 
increased  error  with  increasing  time  step.  This  is  because  the 
truncation  error  for  the  ADI  and  CN  methods  is  of  order  At*  +  h2  ,  and 
the  At2  contribution  is  negligible  for  small  time  steps.  The  IM  and  EX 
methods,  however,  get  contribution  to  the  truncation  error  from  At,  not 
At2  ,  so  it  is  not  negligible.  As  a  result,  at  small  time  steps  the  IM 
and  EX  error  increases  as  time  step  increases.  As  the  time  step 
increases  further,  the  At2  component  of  truncation  error  becomes 
significant,  causing  increasing  error  for  the  ADI  and  CN  solutions  as 
can  be  seen  in  Figures  2-7. 

A  finer  mesh  is  expected  to  yield  a  more  accurate  solution  since  it 
reduces  the  truncation  error.  For  example,  for  the  first  problem  with  a 
time  step  of  At  »  0.01  and  a  mesh  spacing  of  h  =  0.1,  the  order  of  the 
truncation  error  of  the  ADI  method  is  h2  +  At2  *  ( 0 . 1 ) 2  +  ( 0 . 0 1 ) 2  * 
0.0101.  Cutting  the  mesh  size  in  half  gives  error  on  the  order  of 
( 0 . 05) 2  +  (0.01)*  »  0.0026.  The  expected  ratio  of  error,  then,  is  R  = 
0.0101/0.0026  =  3.9.  Figures  2  and  3  show  the  error  ratio  is  actually 
R  =  1.55e-10/3.27e~ll  =  4.7.  Thus  the  actual  improvement  in  error  due 
to  halving  the  mesh  size  is  about  what  was  expected. 

Oscillatory  Solutions 

In  Chapter  3,  stability  was  examined  using  the  coefficient  method. 
A  region  of  stable  oscillations  was  predicted  for  the  ADI  and  CN 
methods.  As  predicted,  the  calculations  resulted  in  ADI  and  CN 
oscillations  for  certain  values  of  time  step.  In  some  cases,  the  errors 
resulting  from  the  oscillations  were  large  enough  to  render  the  solution 
useless  for  most  engineering  applications. 


Problem  #1  Oscillations 


The  first  problem  was  of  limited  usefulness,  but  did  reveal  stable 
oscillations  for  the  CN  and  ADI  solutions.  A  sample  of  these 
oscillations  is  presented  in  Figure  8,  along  with  the  IM  and  analytical 
solutions.  The  decrease  of  error  with  time,  indicating  stability,  is 
documented  in  Figure  9.  At  this  value  of  time  step  and  r,  the  EX  method 
was  unstable.  Also,  these  oscillations  are  not  symmetric  at  all  four 
corners.  All  observed  oscillations  for  this  problem  began  with 
completely  symmetric  oscillations,  but  as  time  progressed,  the  nature  of 
the  problem  (with  the  solution  going  to  zero)  caused  the  unsymmetric 
behavior.  In  all  cases,  the  IM  method  was  stable  and  did  not  oscillate. 

Problem  #2  ADI  Oscillations 


The  second  problem  more  closely  resembles  an  engineering  problem 
and  is  much  more  useful  for  examining  errors  caused  by  oscillations. 
Figure  10  shows  the  ADI  and  analytical  solutions  for  two  different 
times,  with  a  large  value  of  r  for  which  stable  oscillations  are 
predicted.  It  is  clear  that  the  oscillations  decreased  as  time 
progressed. 

Figure  11  confirms  that  the  errors  decrease  with  time,  so  the 
solution  is  shown  to  be  oscillatory  but  stable.  Figure  11  also  shows 
the  oscillations  induced  errors  on  the  order  of  10*  -  100*,  too  great  an 
error  for  most  engineering  applications. 
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Error  Progression  With  Time  Problem  #1 
At=0.1  h-0.1  r=At/hJ=10.0 
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ADI  SOLUTION 


Figure  20.  Numerical  Solutions  Problem  #2 

T ima=0. 1  At=0.0125  h=0.05  r=At/hJ=5.0 
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Figure  18.  Numerical  Solutions  Problem  #2 

Time=0.4  At=0.1  h=0 . 1  r=At/h*=10.0 
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Figure  17.  Error  Progression  With  Time  Problem  #2 
Ata0.025  h-0.1  r=At/h*=2.5 
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ADI  SOLUTION 


CN  SOLUTION 


Figure  16.  Numerical  Solutions  Problem  #2 

T ime*0. 1  a  t=0. 025  h-0.1  r=At/h*  =  2.5 


Figure  16  shows  the  solutions  computed  with  the  ADI,  CN ,  and  IM 
methods  for  a  value  of  r  =  2.5.  Note  that  this  is  in  the  oscillatory 
region,  so  the  unstable  EX  method  gave  no  solution.  As  predicted  for 
this  value  of  r,  the  ADI  and  CN  methods  resulted  in  oscillations  and  the 
IM  method  did  not.  To  see  if  these  oscillations  are  stable,  the  errors 

are  plotted  against  time  in  Figure  17.  Clearly  the  error  decreases  with 
time,  so  all  three  methods  are  stable. 

Using  a  larger  time  step  which  gives  a  value  of  r  *  10,  Figure  18 
again  reveals  oscillatory  behavior  for  the  ADI  and  CN  methods.  Figure 
19  shows  that  the  error  decreases  with  time,  so  the  methods  are  stable. 
For  a  finer  mesh  with  h  *  0.05,  Figure  20  shows  another  example  of 
oscillations,  this  time  with  a  value  of  r  *  5.0.  Figure  21  again 
demonstrates  that  the  methods  are  stable,  as  the  error  decreases  with 
time. 

To  examine  the  effects  of  the  CN  oscillations,  the  maximum  error 
for  the  CN  solutions  for  different  values  of  r  in  the  oscillatory  region 
are  plotted  in  Figures  22  and  23.  As  expected,  the  larger  the  time  step 
(and,  therefore,  the  larger  r  is),  the  larger  the  error.  As  with  the 
ADI  method  (Figures  14  and  15),  these  figures  show  that  in  some  cases 
the  method  can  be  used  despite  oscillations.  Like  the  ADI  method,  if  an 
accuracy  criteria  of  maximum  error  less  than  1*  is  used,  values  of  r 
less  than  two  are  usable,  values  of  r  of  4  -  5  marginal,  and  r  =  10  is 
unusable.  Compared  to  the  ADI  solutions  (Figures  14  and  15),  the  CN  is 
seen  to  have  about  the  same  amount  of  oscillations  as  the  ADI  methods 
for  small  r,  and  greater  error  at  large  r. 
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Figure  14.  ADI  Solutions  Problem  12  h=0.1 
Maximum  Error  vs  Time  For  Various 
Values  Of  r  In  The  Oscillatory  Region 
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Figure  12.  ADI  and  Analytical  Solutions  Problem  #2 
At*0.025  h»0 . 1  r®At/hJ  =  2 .5 


For  a  smaller  time  step  and,  therefore,  a  smaller  value  of  r  (still 
in  the  oscillatory  region).  Figure  12  shows  the  ADI  and  analytical 
solutions  at  two  different  times.  Again,  the  error  decreases  with  time, 
as  can  be  seen  in  Figure  12  and  confirmed  in  Figure  13.  In  this  case, 
the  errors  at  times  greater  than  about  0.3  seconds  are  on  the  order  of 
1*,  and  so  may  be  useful  for  engineering  applications. 

The  ADI  method  also  gives  oscillatory  solutions  for  other  values  of 
r  >  1/2.  The  oscillations  are  sometimes  dominated  by  the  other  errors 
and  cannot  be  observed  in  the  solution.  Generally,  the  greater  the 
value  of  r,  the  more  severe  the  oscillations.  For  all  cases  of 
oscillatory  solutions,  it  was  found  that  the  error  decreased  with  time, 
so  stability  was  maintained. 

Figures  14  and  15  show  the  maximum  error  as  time  progresses  for 
different  values  of  r  for  which  oscillations  are  predicted.  Figure  14 
is  for  a  mesh  spacing  of  h  *  0.1  and  Figure  15  is  for  h  *  0.05.  As  time 
increases,  the  error  decreases  and  the  oscillations  damp  out.  It  is 
seen  in  these  figures,  however,  that  at  early  times  the  error  is  large. 
If  the  requirement  for  accuracy  in  an  engineering  problem  was  for  the 
maximum  error  at  any  point  to  be  less  than  1*,  then  Figures  14  and  15 
show  that  values  of  r  of  2  or  less  may  be  usable,  values  of  r  of4  ■  5 
marginal,  and  r  *  10  unusable.  Also,  the  longer  the  time  of  interest, 
the  larger  the  time  step  (and  thus  larger  r)  that  can  be  used  since  with 
enough  time,  all  the  errors  decrease  substantially. 

CN  and  Oscillatory  Results 

As  discussed  in  Chapter  3,  the  CN  method  may  have  stable 
oscillations  and  the  IM  method  is  expected  to  always  be  stable  without 
oscillations. 
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Figure  11.  Error  Progression  With  Time  ADI  Problem  12 
At-0.1  h-0.1  r*At/hl=10.0 
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Figure  10.  ADI  And  Analytical  Solutions  Problem  12 
At*0.1  h=0.1  r=At/h**10.0 
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Summary  of  Oscillations 

Dozens  of  cases  with  different  mesh  spacing  and  time  steps  were 
run.  In  all  the  cases,  behavior  was  as  predicted.  For  values  of  r 
greater  than  1/2,  both  the  CN  and  ADI  methods  had  stable  oscillations. 
Oscillations  did  not  always  appear,  however;  in  some  cases,  the  dominant 
error  was  truncation  error.  As  r  increased,  so  did  the  degree  of 
oscillations.  In  every  case,  both  oscillatory  and  non-oscill atory ,  the 
CN  and  ADI  solutions  were  stable. 

As  predicted,  the  IM  method  never  oscillated,  even  for  large  values 
of  r.  This  non-oscill atory  behavior  did  not  ensure  accuracy,  however. 
Although  more  accurate  than  the  CN  and  ADI  methods  for  large  time  steps, 
often  the  IM  solution  had  enough  error  to  be  unusable.  In  every  case 
run,  with  values  of  r  up  to  200,  the  IM  method  was  stable.  Of  course, 
since  r  ■  aAt/h2 ,  a  different  value  of  o  will  affect  the  choice  of  time 
step. 


Computer  Time  Required 

Central  processor  unit  (cpu)  time  required  for  calculations  was 
measured.’  For  a  set  number  of  iterations,  the  EX  method,  which  doesn't 
require  solving  a  system  of  simultaneous  equations,  was  the  fastest. 
The  ADI  method,  which  results  in  a  tridiagonal  matrix  system,  was  second 
fastest.  The  finer  the  mesh,  the  greater  the  order  of  the  matrix,  and 
so  the  greater  the  advantage  the  ADI  method  has  over  the  CN  and  IM 
methods.  This  is  seen  in  Figures  24-26,  which  give  the  time  required  to 
solve  representative  problems  with  a  coarse  mesh  spacing  of  h  =  0.1,  a 
finer  mesh  of  h  *  0.05,  and  a  very  fine  mesh  of  h  »  0.025.  The  relative 
time  advantage  of  the  ADI  method  increases  tremendously  as  the  mesh  is 
refined. 
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Figure  26.  CPU  Time  Used  To  Solve  Problem  #2 
20  Time  Steps  of  Ata0.0001 
Mesh  Spacing  h-0.025  Time*0.002 
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The  real  test  of  time,  of  course,  is  the  time  required  to  solve  a 
problem  to  the  desired  level  of  accuracy.  To  test  this,  the  solution  to 
the  second  problem  was  computed  at  0.4  seconds  using  the  largest  time 
step  for  which  the  maximum  error  at  any  grid  point  was  less  than  1*. 
The  computations  were  made  for  mesh  spacings  of  h  =  0.1  and  h  *  0.05. 

For  the  case  of  h  *  0.1,  the  stability  condition  required  160  time 
steps  with  At  *  0.0025  for  the  EX  method.  The  ADI,  CN,  and  IM  methods 
all  reached  the  desired  accuracy  with  16  time  steps  of  At  *  0.025. 
Figure  27  shows  that  the  ADI  method  is  fastest,  and  the  EX  the  slowest. 
For  such  a  course  grid,  the  difference  is  noticeable  but  relatively 
small.  For  a  refined  mesh  spacing  of  h  *  0.05,  the  ADI  proved  far 
superior.  The  matrix  equation  requiring  solution  has  a  much  larger 
order  than  for  h  *  0.1  (361  compared  to  81),  and  so  the  advantage  of  the 
tridiagonal  system  is  greater.  Figure  28  shows  the  cpu  time  required 
for  the  computations.  Oscillatory  errors  limited  the  CN  and  ADI  methods 
to  32  time  steps  with  At  *  0.0125,  and  the  IM  required  only  16  time 
steps  with  At  »  0.025.  Despite  twice  as  many  time  steps,  the  ADI  method 
was  nearly  four  times  faster  than  the  IM  method.  Due  to  the  stability 
condition,  the  EX  method  required  640  time  steps  with  At  =  0.000625, 
and  was  the  slowest  method. 

In  summary,  the  finer  the  mesh,  the  more  speed  advantage  the  ADI 
method  has.  Even  with  limitations  on  the  time  step  due  to  oscillatory 
behavior,  the  ADI  was  the  fastest  method. 

Chiver *s  Thesis 

The  primary  error  in  Chiver's  thesis,  which  was  the  starting  point 
for  this  study,  was  his  misinterpretation  of  the  results.  He  stated 
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that  his  oscillatory  solutions  were  unstable  (3:47),  but  did  not 
demonstrate  that  they  were,  having  failed  to  test  for  stability  by 
studying  error  progression  with  time. 

Chiver's  results  and  conclusions  on  the  accuracy  of  the  methods 
(3:42,47)  are  not  clearly  stated  and  are  not  documented  well  enough  to 
support  his  assertions.  He  also  gives  the  improper  AOI  truncation  error 
(3:22)  with  no  reference. 

Small  errors  were  the  improper  application  of  the  initial  condition 
in  the  derivation  of  the  analytical  solution  (3:5,55)  and  a  missing 
factor  of  two  in  the  Crank-Nicolson  formula  (3:16).  It  is  likely  that 
these  two  errors  were  typographical. 


V.  Conclusions  and  Recommendations 


Conclusions 

The  comparison  of  the  four  finite  difference  methods  did  not  reveal 
any  unexpected  behavior.  It  was  shown  that  certain  stability  and 
oscillatory  requirements  need  to  be  met,  and  the  results  bore  this  out. 
Based  on  the  results  of  this  study,  the  following  conclusions  can  be 
inferred. 

1.  The  ADI  method  will  yield  stable  oscillatory  solutions  with 
larger  time  steps.  The  oscillations  may  create  errors  large  enough  to 
make  the  solution  unusable.  Values  of  r  =  oAt /ax1  *  aAt/Ay*  =  aAt/h* 
greater  than  2  must  be  used  with  caution,  and  values  of  r  greater  than  4 
or  5  create  oscillatory  errors  too  large  to  be  useful. 

2.  The  oscillations  of  the  CN  method  are  similar  to  those  of  the 
ADI  method.  Useful  values  of  r  are  about  the  same  for  the  CN  method  as 
for  the  ADI. 

3.  For  small  time  steps,  the  CN  and  ADI  are  more  accurate  than  the 
IM  method.  The  time  required  to  run  the  CN  and  TM  methods  are  about  the 
same,  so  a  decision  between  the  CN  and  IM  methods  should  consider  the 
time  step  required.  For  large  time  steps,  the  IM  method  was  more 
accurate,  but  the  error  can  still  be  too  large  to  be  useful. 

4.  The  ADI  method,  even  with  restrictions  on  time  step  due  to 
oscillations,  is  faster  than  the  other  methods.  The  finer  the  mesh 
size,  the  greater  the  speed  advantage  of  the  ADI  method. 

5.  As  expected,  the  stability  requirement  severely  limits  the  EX 
method  to  very  small  time  steps. 


Recommendations 


Based  on  experiences  studying  these  methods,  the  following 
recommendations  for  further  study  are  proposed. 

1.  Irregular  geometries  will  permit  the  ADI  to  maintain  its 
tridiagonal  form,  but  the  CN  and  IM  will  not  have  symmetric,  positive 
definite  matrices.  These  methods,  then,  will  have  to  be  solved  with  a 
slower  method  than  Cholesky  (square  root)  decomposition.  A  study  of 
irregular  geometries  may  reveal  an  even  more  significant  speed  advantage 
for  the  ADI  method. 

2.  Two  dimensional  problems  with  Neumann  or  mixed  boundary 
conditions  could  investigate  complications  caused  by  these  boundary 
conditions. 

3.  Comparisons  of  the  ADI  with  Patankar's  exponential  (power) 
method  (9:32)  would  be  useful  for  investigating  this  relatively  new 
method. 

4.  Comparisons  of  the  methods  for  spherical  and  cylindrical 
geometries  may  be  useful,  as  these  are  commonly  encountered  geometries 
in  engineering. 

5.  A  comparison  of  the  ADI  with  other  methods  in  three  dimensions 
and  a  study  of  oscillatory  behavior  for  three  dimensions  could  give 
further  insight  into  the  ADI  method. 

6.  The  stability  conditions  (analyzed  here  using  the  coefficient 
method)  should  be  extended  by  using  the  matrix  approach  such  as  Myers 
17:284-286)  presents  for  two  dimensions. 
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